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Abstract. We theoretically study the adiabatic preparation of an antiferro- 
magnetic phase in a mixed Mott insulator of two bosonic atom species in a one- 
dimensional optical lattice. In such a system one can engineer a tunable parabolic 
inhomogeneity by controlling the difference of the trapping potentials felt by the 
two species. Using numerical simulations we predict that a finite parabolic po- 
tential can assist the adiabatic preparation of the antiferromagnet. The optimal 
strength of the parabolic inhomogeneity depends sensitively on the number im- 
balance between the two species. We also find that during the preparation finite 
size effects will play a crucial role for a system of realistic size. The experiment 
that we propose can be realized, for example, using atomic mixtures of Rubidium 
87 with Potassium 41 or Ytterbium 168 with Ytterbium 174. 
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1. Introduction 

The hardware of an (analog) quantum simulator 1\ is a controlled quantum system 
that is a clean and tunable realization of a many-body model system of interest 
(see also Refs. |3]). In a quantum simulation such a quantum machine is used 
to experimentally measure dynamical or equilibrium properties of the model that are 
hard to obtain by using a classical machine. A typical protocol for studying equilibrium 
properties will start from a parameter regime of the model that is well understood 
theoretically and, thus, allows validating that a state close to thermal equilibrium 
can be prepared faithfully. In a next step, the system is then guided slowly into the 
parameter regime of interest. If it can be assumed that the dynamics during this 
parameter variation is close to adiabatic, the system will finally be in a state close 
to the target state, which is defined to be the thermal equilibrium for the new set of 
parameters characterized, e.g., by the same entropy and particle number as the initial 
state. 

However, typically a phase transition is expected to occur on the way between 
the initial and the target regime. Crossing this transition is potentially a source for 
an increased production of excitations as described by the Kibble-Zurek mechanism 
in the case of a continuous phase transition (see Refs. [11 [5] and References therein). 
In order to keep defect creation at a minimal level, it has been proposed to bring the 
system from one quantum phase to the other without passing a true phase transition 
by employing spatial inhomogeneity [51 [5]. Namely, in an inhomogeneous system the 
transition from one phase to another can happen as a crossover at a spatial boundary 
of finite width. By parameter variation this boundary can be moved through the 
system at a finite speed, eventually bringing it from one phase to the other. Such 
a strategy is similar to methods like growing crystals or pulling them out of the 
melt. It has been investigated theoretically in the simple model system of a quantum 
spin-1/2 chain (with Ising or XY coupling) in an inhomogeneous transverse field; 
here the ferromagnet-to-paramagnet transition (a continuous phase transition in the 
uniform system) can be induced practically without defect creation, provided the phase 
boundary moves slow enough [7J [H] . 

In this paper we investigate theoretically a problem of direct experimental 
relevance. Namely whether a parabolic inhomogeneity can be useful to assist the 
adiabatic preparation of an antiferromagnetic quantum phase in an experiment with 
ultra cold atoms [9l 110) . The antiferromagnetic crystalline order shall be grown in 
space from the center of the system outwards. To this end, we consider a two-species 
mixture of ultra cold bosonic atoms in an optical lattice with strong on-site repulsion 
(see, e.g., Refs. [TTJ QH US HI E2 EEE1 H3 US] for the equilibirum properties of such a 
system). The system can be described by a quantum XXZ-spin-1/2 model [TTlfT2l[T^] 
and we are interested in the transition from an easy-plane ferromagnetic phase in 
the spin xy plane to an easy-axis antiferromagnetic phase in z direction. We will 
concentrate on a one-dimensional system with the dynamics along the perpendicular 
directions frozen out by a strong transversal confinement. 

The motivation for the present work is twofold. First of all, the quantum 
antiferromagnetically ordered target state is known to be very fragile with respect 
to thermal fluctuations, because it is stabilized by low-energy superexchange physics 
only |17] (for antiferromagnetism in ultra cold atoms without superexchange cf. Refs. 
P33 [23 1211 HI])- This makes its experimental realization challenging. It is, thus, 
desirable and of immediate relevance for current experimental studies to investigate 
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how the state can be prepared with a minimum of heating. A related problem 
of great importance is the preparation of the fermionic Heisenberg antiferromagnet 
being a prerequisite for mimicking the intriguing physics of high-temperature cuprate 
superconductors [23] with ultra cold atoms [23] . 

Another motivation lies in the fact that the system we are studying possesses 
several interesting properties. It allows to experimentally control and (despite of the 
fact that the particles are always trapped) even switch off completely a parabolic 
inhomogeneity by tuning the relative trap strength of the two bosonic species [25] . 
This enables the experimentalist to study the influence of (in)homogeneity in detail 
also in the laboratory. 

The system is also rich and generic enough to give rise to effects that potentially 
disfavor the usage of inhomogeneity for the purpose of an adiabatic state preparation. 
For example, mass flow can be a limiting factor, especially if domains of insulating 
phases appear, acting as barriers that hamper density redistribution, as recently 
discussed in the context of the inhomogeneous bosonic Mott transition [25] [27] , 
Another aspect is that the transition can change from a continuous (second order) 
transition to a discontinuous (first order) transition (see Ref. 28 for the two- 
dimensional case). The continuous transition occurs without inhomogeneity, whereas 
the discontinuous transition is relevant for studying the transition with inhomogeneity. 

For the related problem of a fermionic Heisenberg antiferromagnet adiabatic 
protocols based on inhomogeneities that reduce the discrete translational symmetry 
of the system, have been investigated recently [29 ] [30 ] , 

In this paper, we study the influence of a static parabolic inhomogeneity, while 
the transition from one quantum phase to the other is induced by varying terms that 
themselves do not break translational symmetry. Similar scenarios have recently been 
investigated for temperature-driven phase transitions [3T] [32J [33] or for ramps not 
passing a phase transition [34] . We are not considering a quench of the inhomogeneity 
itself, as it has been studied for example in Ref. [35] . We note in passing that 
the dynamics of a harmonically trapped bose-bose mixture in response to a sudden 
displacement of the trap has recently been investigated theoretically as a probe for 
different quantum phases in the system |3"S] . 

The numerical studies of the one-dimensional system that we present in this paper 
indicate that a parabolic potential can indeed assist the adiabatic preparation of the 
antiferromagnetic target state. However, the optimal strength of the inhomogeneity 
depends in a sensitive way on the imbalance between the two bosonic species; the larger 
the imbalance the larger the optimal inhomogeneity. Therefore, the possibility to tune 
the inhomogeneity [25J should be an advantage for preparing the antiferromagnetic 
order. We also observe that for realistic system sizes the time evolution during the 
parameter ramp into the antiferromagnetic regime is still governed by finite size effects 
that go beyond the local-density picture. Namely we find precursing antiferromagnetic 
correlations already in the ground state of the system outside the antiferromagnetic 
regime. These are contaminated with imperfections (like kinks) that originate from the 
inhomogeneity. The imperfect correlations are amplified when the system is ramped 
into the antiferromagnetic regime. It is difficult for the system to get rid of these 
imperfections, as it would be required for a perfectly adiabatic time evolution. So we 
find the best results for parameters giving rise to an initial state with a low degree of 
imperfections in the precursing antiferromagnetic order. 

The paper is organized as follows. The system and the different models describing 
it are introduced in section [5] The structure of the grand canonical ground state 
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phase diagram is reviewed in section [3] with some details on how the phase diagram 
has been computed using the Bethe ansatz in appendix | Appendix A| The protocol for 
the preparation of the antiferromagnetic state is described in detail and motivated in 
section[4j The results of our numerical simulation of this protocol are finally presented 
in section [5] before we conclude in section [6] 



2. System and models 

We are considering a system of ultra cold atoms given by a mixture of two bosonic 
species in one spatial dimension (ID) subjected to a steep optical lattice potential. In 
recent experiments such mixtures have been loaded into optical lattices, among them 
Potassium (K41) Rubidium (Rb87) mixtures [37] and mixtures of different hypcrfinc 
("spin") states of Rb87 [33] ESI SOI S3 S3 S3 ■ Other candidates include mixtures of 
different Ytterbium-Isotopes [H] SS] that offer a rich variety of scattering properties 
depending on the selection of isotopes [IB] . In the following we consider a two-species 
system characterized by all-repulsive interactions and with the intraspecies repulsion 
being strong compared to the interspecies repulsion. Such a system can be realized 
experimentally by using an Ybl68-Ybl74 mixture [46j or, alternatively, by taking a 
K41-Rb87 mixture with the interspecies scattering length tuned small by means of a 
Feshbach-resonance [27]. 

The ID mixture of two bosonic species s — a, b in an optical lattice is described 
by the Bose-Hubbard Hamiltonian 



oo 

#bh = J2 J s{a\ l+1 a sl 



+ (V s i - Ms)^} + U ab h at h u 



+ 1 ^h a i(h al > - 1) + -^-hu{nu - 1) i (1) 

where a se and n s e are the bosonic annihilation and number operator for particles of 
species s at lattice site £. Tunneling between neighboring sites is captured by the 
positive matrix elements J s ; the three positive Hubbard energies U a b, U aa , and Ubb 
characterize the repulsive inter and intra species on-site interactions; and the particles 
are confined by the harmonic trapping potentials V s e — ^a s £ 2 . In the ground state the 
total numbers N a and Nf, of a and b particles are controlled by the chemical potentials 

Ms- 

We are interested in the regime of strong repulsive interactions with the Hubbard 
energies U s / S being positive and large compared to both the tunneling matrix elements 
J s and the chemical potentials fi s such that double occupancy is strongly suppressed. 
Under these conditions we can effectively describe the system within the low-energy 
subspace S defined by 

S : n a i + h M < 1 W. (2) 

We employ degenerate-state perturbation theory [48] up to second order with respect 
to tunneling processes and obtain the effective Hamiltonian 



^(V;^ - fi s )n se - J s P s (al i+1 a si + h.c.)P s 

(3) 



J a s i + x a M+x a M a st ~ ^n se+1 n se 
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acting in S. Here the operator P$ projects into the subspace S and s denotes the 
species opposite to s. The first and second term originate from zeroth- and first-order 
perturbation theory, respectively. The new matrix elements J and U describe second- 
order superexchange processes. While J quantifies swaps between a and b particles on 
neighboring sites, U characterizes an attractive nearest-neighbor interaction between 
a and b particles. These matrix elements stem from perturbative admixtures of Fock 
states with one site occupied by both an a and a b particle and read 

J = 2^ (4) 

Uab 

U = 2±f±. (5) 

Uab 

Furthermore, when deriving H e g two additional approximations have been made 
for simplicity that both are well justified. The first one is that we neglected second- 
order terms involving virtual excitations (perturbative admixtures) with two particles 
of the same species on the same site. The amplitudes of such terms are proportional 
t° Ja(b)/U a a(bb) an( I are much smaller than the effective matrix elements fit) and (51), 
since we assume U a b <C U aa , Ubb- The second simplification consists in neglecting ffie 
small potential energy differences between neighboring sites (V^+i— V s e) = a s (£+l/2) 
that would appear together with U a b m the denominators of the second-order matrix 
elements Q and ([5]). This approximation is well justified for typical slowly varying 
traps. 

Concerning the level of approximation, the description of the bosonic system in 
terms of the effective Hamiltonian ^ is comparable to the iJ-model for spin-f/2 
fermions on a lattice with strong on-site repulsion [33 HO] • It describes the low energy 
physics of a doped magnet by combing two elements; the superexchange coupling 
between the spin (or species) degree of freedom on neighboring occupied sites on the 
one hand and, on the other hand, the dynamics of the charge (or total density) degree 
of freedom due to the presence of holes (vacant sites). For fermions the interplay 
between both is conjectured to give rise to intriguing physics like high-temperature 
superconductivity in the case of a square lattice (23] . The homogeneous version of the 
bosonic model ^ has been studied theoretically in Refs. [5TJE21E3] where, e.g., phase 
separation between hole-rich and hole-free regions is predicted on the square lattice. 

For slowly varying traps it is useful to introduce the local chemical potentials 
Hsl = fi s — V s i- For sufficiently large \i a and [Xb (i.e. for a sufficiently large total 
particle number N a +Nb), in an extended region M in the trap center the local chemical 
potentials will be large enough to strongly suppress the existence of unoccupied sites. 
(An estimate of the size of M will be given at the end of this section) . In this region the 
particles form a mixed Mott insulator with occupation (n a + hb) — 1. The remaining 
degrees of freedom, namely which site is occupied by which species, can then effectively 
be described within the subspace S' of unit filling, 

S' : n ae + h e = 1 V£e M. (6) 

In 5" and for I € M the Hamiltonian is again given by H c s, but the tunneling terms 
can now be dropped, giving 

H = Yj [ ~ J ^ai+i h M+Ai h ai + h - c -) 
t 

- -r("M+l"«< + n al+lflu) + - [Vi - H)(hat ~ flbl) 
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We have also omitted the constant terms \{^ a l + MmX^W + n w) — \{li*al + Mm) an< i 
introduced the notation 

V t = Vai ~Vu = \ («« - a b )£ 2 = X -at 2 (8) 

and 

= fx a - /j, b . (9) 

The effective Hamiltonian Q will be the starting point for the remaining sections of 
this paper. 

The inhomogeneity V# appearing in H is characterized by the difference of the 
trap frequencies a — a a — ab- This can be explained as a consequence of the constraint 
h a i + hbi = 1; tunneling of an a particle from t to I' has to be combined with the 
counterflow of a b particle from i! to i. In an experiment the degree of inhomogeneity 
a can be tuned continuously, simply by adjusting the trapping potentials of the two 
species with respect to each other. In particular, the accessible parameter space 
contains the homogeneous model with a = that is realized for equal traps a± = oti 
[2"5] (as well as the regime of a < 0) . The fact that this limit can be reached without 
the model description breaking down is a crucial ingredient of the experiment proposed 
here. 

Apart from the dimensionless inhomogeneity a/U, the Hamiltonian H describing 
the simulator region M of our system is characterized by two further dimensionless 
parameters, namely J/U = (J a /Jb + Jb/Ja) 1 and jJb/U. In an experiment these 
can be controlled independently by adjusting the ratio of tunneling strengths J a /Jb 
(controlled by the lattice depths for a and b particles) and the imbalance between a 
and b particles (N a - N b )/(N a + N b ). 

It is both instructive and convenient to express the Hamiltonian ([7]) in terms 
of composite-particle [54] and spin (13j degrees of freedom. The former description 
is obtained by introducing composite particles that are hard-core bosons with 
annihilation operators 

K = ^bl^at ( 10 ) 

and number operators fi£ = SJS», such that [b», b\,] — #^'(1 — 2n^) and 

hk = = (11) 

in S' . With that ne = b\b t < 1. In S' the composite-particle occupation numbers are 
equal to those of the a particles, hi = h a i(n b i + 1) = n a , whereas "composite holes" 
correspond to b particles, 1 — hi = h b . The Hamiltonian Q can now be re- written 
like 

H 



E 



[ - J(bl + A + h.c.) + U (ne+i - -) (hi - -) 
+ (y t -n)(ht-±)]. (12) 

This Hamiltonian describes hard-core bosons in a tunable trapping potential Ve , with 
hopping matrix element J, repulsive nearest neighbor interactions U, and chemical 
potential \x. 

A spin- 1/2 description is defined by identifying the species s with an internal spin 
degree of freedom with spin f (I) for a (b) particles. Introducing the vector of Pauli 



Quantum crystal growing 



7 



matrices cr s > s for this spin degree of freedom, we can define the spin operator at site 
I: 



S,= 



\Yl^sU cr e'sa si (13) 



2 

S ! S 



with components 

Sf = l(al e a M + h.c.) (14) 

% = l&h& ai +h.c) (15) 

S$ = -(fi a e - flu)- (16) 

In the subspace S' these spin operators Se describe a spin- 1/2 degree of freedom at 
every site. In terms of these degrees of freedom the Hamiltonian takes the form 

£ = £[- 2J0s?s? +1 + sjs* +1 ) + usfst + <y t n)si] (17) 
i 

of an XXZ spin chain with ferromagnetic spin coupling —2J = J x = J y in x 
and y direction, antiferromagnetic Ising coupling +U = J z in z direction, and an 
inhomogeneous magnetic field (Vi — fx) = hi in 2 direction. 

In the following we will focus on the central Mott insulator region M described by 
the Hamiltonian H. It serves as a simulator for the dynamics described by the model 
Hamiltonian H with tunable parabolic inhomogeneity. Let us briefly estimate the size 
of the Mott insulator region in the zero-temperature equilibrium state. In the limit of 
zero tunneling both species fill up the trap such that sites I with \£\ < (N a +Nb)/2 = £ 
are occupied. If fi s * is the larger of the two chemical potentials \x a and [if,, and a s * 
denotes the corresponding trap frequency, one has |a s *£g = fi s * . In order to suppress 
doubly occupied sites near the trap center besides J s <C V a b also fi s * < U a b is required. 
The latter implies that 2£ < ^8U a b/a s *. 

When finite tunneling is included, the edge of the occupied region will soften. 
With increasing \£\, near \£\ = £ the occupation (n a + hb) will drop from 1 to 
within a crossover region of width A£. The width A£ is such that the increase of 
the trapping potential in the crossover region is of the order of the tunneling matrix 
element. More precisely, diV s t\i=t a A£ s — ae£oA£ s ~ J s gives A£ s ~ J s /(a s £o) and 
one has either A£ — max(A^ a , A£b) if both species can be found near the edge \£\ w £q 
(that is \i a ~ Mb) or ^ — if basically only s* particles occupy the edge region 
(i.e. [A s » <C Ms*)- l n the former (more restrictive) case the crossover region A£ will be 
small compared to £q as long as J s <C ot s £\ ~ 2/z s * < 2U a b, which has been required 
already. Therefore, the "simulator region" M has an extent of L ~ 2£q — 2A£ that 
will be of the order of £q and it can host an extensive fraction of the particles in the 
system. 



3. Ground state phase diagram of the homogeneous system 



In the central Mott region M the system is described by the Hamiltonian H that we 
expressed in different representations. In the following we will stick to the language of 
the hard-core boson model (12), unless we explicitly mention the two-species [Eq. Q] 
or spin [Eq. ( 17)] description. For a homogeneous system with Ve. — 0, the ground state 
of this model is characterized by two dimensionless parameters, the scaled chemical 
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Figure 1. Ground-state phase diagram of hard-core bosons in a one-di mens ional 
lattice with nearest-neighbor repulsion U described by the Hamiltonian i ] 1 2[ ) with 
V( = 0. Subfigure (b) shows the parameter plane given by the dimensionless 
chemical potential fj,/U and tunneling parameter J/U. Thick solid (dashed) lines 
indicate continuous (discontinuous) phase transitions. The system possesses a 
compressible superfluid phase (SF) and three incompressible phases, the vacuum 
at zero filling (V), the particle-hole reflected vacuum at unit filling (U), and a 
density- wave insulator at half filling (DW). The color code describes the filling 
n and clearly indicates the incompressible regions where dn/dfi = 0. The DW 
insulator is characterized by a long-ranged staggered order in the site occupation, 
measured by the order parameter Odw defined in Eq. jl8[ l . Odw is non-zero only 
in the DW phase, where it depends on J/U, as plotted in subfigure (a), but not 
on /i/U. The sudden jump of Odw to zero when entering the SF phase renders 
the DW-to-SF transition discontinuous (though the filling n varies continuously 
when passing the transition). The only exception occurs at the tip of the lobe- 
shaped DW domain (located at fi/U = and J/U = 1/2) where the transition is 
continuous. The different thin curves plotted in (b) indicate the chemical potential 
/i/U of finite systems in the local density approximation with N = 31,29,27,25 
particles on L = 61 sites subjected to a trap of strength a/U = 10 -3 . The vertical 
"system lines" attached to these curves indicate how much the chemical potential 
drops between the center and the edge of the system according to the local density 
approximation. 
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potential fi/U and the scaled tunneling matrix element J/U, with the nearest-neighbor 
repulsion U serving as the unit of energy. We have computed the phase diagram 
of the homogeneous system in the fi/U-J/U plane by employing the Bethe-Ansatz 
solution developed in Refs. [55J [53 [57] . Details of this solution are given in appendix 
Appendix A In Fig. [ljb) we plot the zero-temperature phase diagram and the boson 



filling n = j- 2^2 eeM (hi) with L = X^eA/ ^ denoting the number of lattice sites in M. 

The phase diagram shown in Fig. |T|(b) possesses the following structure. First of 
all, it reflects the particle-hole symmetry of the homogeneous hard-core boson model 
|l2| ; replacing b e — > b\ [implying (hi — 1/2) — > (1/2 — hi)] and fi — > —fi leaves the 
Hamiltonian unchanged, such that fi — > —fi simply interchanges the role of particles 
and holes. 

The energy of a single particle with respect to the vacuum energy is given by 
—fj, — U — 2 J with the kinetic energy reduction — 2J stemming from derealization. 
Therefore, below a chemical potential of fi v /U — —2J/U — 1 the system is in the 
vacuum (V) state \v) with no particles present. Accordingly, for chemical potentials 
larger than fi u /U — —fi v /U — 2J/U + 1 the ground state is the particle-hole reflected 
vacuum, that is the incompressible insulating state \u) = n^j?l u ) & t umt filling (U) 
with exactly one hard-core particle at every site. 

Starting from the vacuum state and increasing the chemical potential fi/U, for 
non-zero tunneling J/U the particle number starts to grow in a continuous fashion 
once the critical parameter fi v /U is passed. Here the system enters a superfluid 
(SF) phase in a second-order phase transition. This phase is characterized by a finite 
compressibility d^n ^ 0, a homogeneous density distribution n£ = (hi) = n, and quasi- 
long-range off-diagonal order, i.e. the correlation function (b\b t ) decays algebraically 
for large \£-£'\. 

For both chemical potential fi/U and tunneling J/U small, another incompressible 
phase is found at half filling, a density wave (DW) Mott insulator [see Fig.jljb)]. This 
phase can be understood by starting from the limit of zero tunneling J/U = 0. Here 
for — 1 < fi/U < 1 a DW state with one particle on every other site is favored as ground 
state. This state is two-fold degenerate and breaks the translational symmetry of the 
lattice. It possesses an energy gap A = min(A p , A^) where A p = U —fi and A^ = fi—U 
are the energy costs for adding a particle or adding a hole (removing a particle) 
somewhere in the system, respectively. [Particle number conserving particle-hole 
excitations (created e.g. if one particle tunnels to a neighbori ng site) come with the 



larger energy cost A = A p + A/j = 2U.] Since the Hamiltonian ( 12 ) conserves the total 
particle number the gap A protects a state at half filling from competing states with 
different particle numbers also for finite tunneling J/U , roughly as long as A is larger 
than the derealization energy 2 J. The rough estimate A = 2 J for the phase boundary 
(corresponding to first order perturbation theory with respect to tunneling) explains 
the lobe shape of the DW insulator domain in the phase diagram and (accidentally) 
even gives the correct critical tunneling strength (J/U) c = 1/2 at the tip of the lobe. 
Actually, this value of 1/2 is fixed by symmetry, namely as the Heisenberg point of 
the model in spin representation [Eq. (17)] where \J X \ = \J y \ — \ J Z \M 



| Re-defining Si — > — Si on every other site also the Ising coupling in z-direction becomes 
ferromagnetic. Now the DW order corresponds to (long-range) ferromagnetic order in z direction 
and the SF phase to (quasi-long-range) ferromagentic order in the xy plane. At the Heisenberg 
point the system is known to possess spin-isotropic quasi-long-range ferromagnetic order. Now 
increasing/decreasing the z-coupling relative to the xj/-coupling slightly, an easy axis/plane is created 
that immediately attracts (at least part of) the ferromagnetic correlations, guaranteeing DW/SF 
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Within the DW domain the particle number and with that the whole structure of 
the ground state does not depend on the chemical potential \i/U . Thus, (within this 
domain) also the DW order is a function of J/U only. It can be quantified in terms 
of the long-ranged density-density correlations by using the order parameter 

Odw= lim ((2n e -l)(2h e , -l))(-iy~ e ' (18) 

\£ — £' | — yoo 

that assumes values between and 1. An analytical expression [58j§]is given by 

oo 

Odw = tanh 4 (nA), cosh(A) = U/2J (19) 

n=l 

and plotted in Fig. [l}a). The fact that the order parameter Odw depends on J/U 
only implies immediately that Odw drops in a discontinuous fashion from a finite 
value to zero when the phase boundary of the DW domain is passed. This makes 
the transition from DW to SF first order almost everywhere. As the only exception, 
the DW-SF transition becomes continuous (second order) through the tip of the lobe; 
this describes the transition at fixed particle number (half filling) driven by tunneling. 
Note that unlike the case of a two-dimensional square lattice with true long-range 
order in the SF phase |28j . where also the particle number varies in a discontinuous 
fashion at the DW-SF transition, in one dimension the filling (n) continuously departs 
from 1/2 when entering the SF phase. 

All in all, at T — the model possesses four different phases. A homogeneous SF 
phase and three distinct insulating phases, the vacuum (V) with ni — 0, the particle- 
hole reflected vacuum with unity filling (U) ni = 1, and the DW Mott insulator at 
half filling n = 1/2 with alternating site occupations. In the original two-species 
picture Q, the n = 1 and n = insulators correspond to a Mott insulator state 
consisting solely of a or 6 particles, respectively; the DW insulator is a Mott insulator 
with a staggered configuration of a and b particles, and the SF phase is a counterflow 
superfluid where a superflow of a particles is accompanied by the corresponding back 



flow of b particles such that (h a +hb) = 1. Finally, in the spin language (17) the n = 
and n = 1 insulators correspond to the fully z-polarized state, the DW state is a phase 
with antiferromagnetic long-range order in the z-components of the spins, and the SF 
state corresponds to quasi-long-range ferromagnetic order in the ccy-plane. 

4. Protocol: Quantum Crystal Growing 

Starting in the SF regime, we wish to study the adiabatic preparation of the crystal- 
like DW insulator state by slowly lowering the tunneling parameter J/U. In particular, 
we are interested in the role played by an inhomogeneity in the form of a parabolic 
potential Vi during this process. For this purpose we consider a finite system of 



L sites described by the Hamiltonian ( 12 ) and characterized by the number of 



hard-core bosons N and by the scaled trap frequency a/U. We mimic the finite 
extent of the simulator region by employing open boundary conditions such that 
I = —R,—R + 1,...,R with R — (L — l)/2. Initially, the tunneling parameter 
J/U assumes a finite value (J/U)q and the system is prepared in its ground state. 

order. (This argument does not exclude an intermediate supersolid phase with both orders present, 
however such a phase is not found within the Bethe ansatz solution.) 

§ There is a misprint in |58l p. 186]. In the formula immediately preceding Eq. (245) a 2 should be 
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Then J/U is ramped down to zero at constant rate within a time span of duration 
T = tH/U. In order to quantify the degree of adiabaticity, after this ramp the degree 
of DW order is measured M 

In order to motivate such a protocol and to gain intuition for the physics related 
to the presence of the parabolic potential, it is instructive to discuss the protocol 
described in the preceding paragraph in terms of the local density approximation 
(LDA). Introducing the local chemical potential m = \i — Vi one assumes that 
the ground state of the inhomogeneous system can locally be approximated by the 
properties of the homogeneous system (summarized in the phase diagram of Fig. [T]) 
with the chemical potential given by /^ffl Within the picture of the LDA, the state of 
an inhomogeneous system with tunneling J/U is represented by a vertical line of finite 
length (the "system line" ) that cuts through the phase diagram of Fig. [I] One end of 
this line lies at fig=o/U — fi/U and corresponds to the center of the trap. The other 
end, to be identified with the edges of the system, lies at (j,£=± R /U = fi/U—aR 2 /(2U). 
So the length of the system line A/i = (a/U)R 2 /2 is directly proportional to a/U. In 
the following we will always assume that a > such that the upper end of the system 
line corresponds to the trap center; results for a < can be inferred from particle-hole 
reflection. 

The chemical potential fi is determined such that the total number of particles 
in the system is given by N. That is, the system line simply shifts upwards when 
the particle number is increased. When J/U is varied, the chemical potential \x has 
to be adjusted in order to keep the particle number N fixed. So when we think 
of adiabatically decreasing J/U the system line will move not only leftwards, but is 
displaced also in vertical direction. In Fig. |T|(b) this is exemplified for three different 
sets of parameters. The non-solid lines indicate how fi/U changes with J/U when the 
particle number is fixed. The vertical lines attached to these lines indicate the system 
line. 

There are good reasons to expect that the presence of a parabolic inhomogeneity 
might be helpful for the adiabatic preparation of the target state (the DW crystal at 
J/U = 0). Consider a slow parameter variation following the protocol that is described 
by the short-dashed thin line in Fig. [ijb). When J/U is lowered the transition to 
the DW phase happens first at the center of the trap (corresponding to the upper 
end of the system line), roughly near J/U — 0.15. From then on, the symmetry- 
broken DW structure can smoothly grow from the center outwards. This process 
resembles the physics of growing a crystal or pulling it out of the melt. However, 
here, crystallization is not driven thermally by lowering the temperature, but rather 
by quantum fluctuations when ramping down the tunneling J/U. Hence, one might 
dub this scheme quantum crystal growing. 

Growing the DW phase in the inhomogeneous system in this way does not involve 
a sharp phase transition (cf. Ref. [3] and references therein). Beyond the local density 
approximation the DW state has a smooth boundary in space. When J/U is lowered, 

| Before evaluating the state and after ramping down the tunneling, one might want to add a further 
step to this protocol in which a/U is ramped down to a/U = 0, such that eventually the system 
becomes homogeneous for all protocols. However, such a step can be omitted; it is irrelevant since at 
J/U = it will not alter the DW order anymore. 

1 One condition for the LDA to be valid is that the variation of the trapping potential from site to 
site should be small compared to the tunneling matrix element J, such that particles can delocalizc 
over larger distances (ten sites, say). For J/U ~ 1 this leads to the requirement (a/U)R 1. On 
the other hand, the healing length, the length scale on which a local perturbation influences the 
many-body wave-function, should be short compared to the spatial structure of the potential. 
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this boundary continuously moves through the system such that the symmetry broken 
crystalline order can growQ 

In the presence of the parabolic inhomogeneity the transition is stretched over 
a finite interval both in the parameter J/U and in time. So neither an accurate 
experimental parameter control nor a precise knowledge of the critical parameter are 
required to control this process. In contrast, for sufficiently large homogeneous system 
the transition happens rather suddenly during the ramp when the phase boundary is 
passed (and it can be expected that the symmetry breaking happens independently 
in remote places of the system such that defects are created). 

Another advantage of the presence of a parabolic potential is that it allows to 
form a crystal in the center of the trap also for particle numbers below half filling. 
The extent Ldw of the DW crystal will depend on the filling n and can be smaller 
than the extent of the full system, Ldw = 2nL < L. In contrast a uniform system 
away from half filling does not possess a DW phase. 

However, we can also anticipate effects that are not in favor of making the system 
inhomogeneous. For example a finite trap (a > 0) necessarily requires filling below 
1/2, which limits the extent of the DW crystal to values Ldw < L. There are two 
different mechanisms that lead to such a constraint. The first one is connected to 
the fact that in the superfluid regime the denisty in the center of the trap will be 
larger than the average density N/L. So when J/U is ramped down, for n = 1/2 the 
DW order will not emerge in the trap center but rather independently at those two 
points (left and right from the center) where the local filling is given by 1/2. As a 
consequence practically no correlation between the crystalline order in both sides of 
the trap will be established (a further detrimental effect connected to this scenario will 
be discussed below) . Filling factors Nj L that avoid this unwanted effect will be lower 
than 1/2 and such that at J/U = 1/2 the local density stays below 1/2 everywhere 
in the trap [this is roughly the case for the protocols with N/L < 0.48 depicted in 
Fig. 0b)]. 

The second mechanism limiting Low * s that the ground state at J/U = at half 
filling can only be a pure DW crystal of length Ldw = L, if the overall potential energy 
drop within the simulator region, Afi — aL 2 /8 (the length of the system line), stays 
below 2U (the width of the DW lobe at J/U = 0, see Fig. [IJ. For larger potential 
drops at the edges and in the center of the trap, regions of vacuum or unit filling 
will form, respectively. In order to avoid a core of unit filling in the center of the 
trap also when (a/U)L 2 > 16, the particle number has to be reduced such that that 
L DW = 2N < 4/yfcT/U. 

A limitation for achieving an almost adiabatic time evolution in the presence of 
a trap can also be given by mass transport. Whereas for the initial state the density 
decreases smoothly from the center of the trap to the edge, the target state possesses 
a DW plateau with a filling of 1/2 (i.e. with one particle per pair of neighboring sites) 
in the center for \£\ < Ldw/2 and a filling of zero for \£\ > Ldw/2. Thus, in order 
to reach the target state the particle density has to be redistributed. Therefore a 
new time scale enters when inhomogeneity is introduced to the system that is not 
related to the physics of the phase transition, namely the time needed to achieve this 
redistribution. 

This is strikingly evident when the mass flow required in order to achieve the 

+ For simple model systems it was found that a sufficiently slow parameter variation guarantees an 
almost almost adiabatic time evolution in such a scenario [7] [8] 
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target state is strongly suppressed by the formation of an insulating domain. The 
protocol at half filling (dash-dotted line in Fig. [TJ is an example for such a situation. 
When insulating DW domains form and grow at two points in the trap, these domains 
divide the system into an inner and two outer regions and they become barriers for 
mass flow between these regions. This means that an adiabatic preparation of the 
target state would require an extremely long ramp time. This detrimental mechanism 
has recently been investigated in the context of the Bose-Hubbard model [26l [27] . 
Note that also when no insulating barrier appears mass flow can still be a factor that 
determines the time required for the adiabatic preparation of the target state. 

In the uniform system (a/U = and the system line shrinks to a point) at half 
filling the transition from the SF to the DW phase happens at the tip of the DW lobe 
and is of second order. For a finite harmonic potential, in turn, according to the LDA 
most parts of the system enter the DW phase at a local chemical potential fig ^ 
and therefore at a tunneling parameter J/U < 1/2 for which the transition is of first 
order in the (grand canonical) uniform system. Of course, corrections to the LDA, as 
we discussed them already, guarantee that in the presence of the harmonic potential 
the phase transition is smoothened into a crossover in space (and also in time when 
the spatial crossover region moves). Nevertheless the crossover will be determined 
by the nature of the phase transition it stems from. We can expect that the larger 
the discontinuity of the first order transition in the uniform system [quantified by 
the jump of the order parameter Odw plotted in Fig. [IJa)] the sharper will be the 
spatial crossover and the smaller will be the rate at which J/U can be changed without 
significantly exciting the system. With respect to this effect, steep traps and low filling 
N/L are not advantageous for an adiabatic preparation of the target state. 

5. Simulation of the time evolution 

In the preceding section we have identified and discussed different mechanisms that 
might play a role when slowly ramping the system from the SF into the DW regime in 
the presence of a parabolic inhomogeneity. While some of them favor the presence of 
the inhomogeneity for the preparation of the DW crystal, others disfavor it. In order 
to find out whether (or when) inhomogeneity has a positive or negative influence with 
respect to adiabaticity, we have simulated the protocol described above numerically 
by using the time-dependent matrix product state ansatz |59[ 160] . 

We consider a realistic system with an odd number of particles N ranging from 17 
to 31 on L = 61 sites with open boundary conditions. These odd numbers are of course 
not crucial, but they ensure that the degeneracy between different symmetry broken 
DW patterns is slightly lifted by the parabolic potential such that our simulation 
always leads to a unique reflection symmetric pattern with larger site occupation at 
the even sites I — 0, ±2, ±4, . . . , ±(JV — 1). Moreover, also in the absence of the 
parabolic potential an odd number of sites L guarantees a unique non-degenerate DW 
ground state at "half filling" N = (L + l)/2 (such that the DW pattern increases the 
occupation on both edge sites). 

In our simulations we compare results for parabolic potentials of four different 
strengths, a/U = 1Q~ 4 , 10 _3 ,6 • 10 _3 ,1Q _2 . We do not switch off the parabolic 
potential completely, since keeping a small finite potential is required for having 
a DW phase also away from half filling. For the two largest values of a/U (the 
two steepest potentials), we only consider particle numbers N of up to 25 and 19, 
respectively, that are smaller than 2/^/a/U. This guarantees that the ground state 
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a/U 


Maximum N 


Afi/U 


aR/U 




1 x 10" 


-4 


31 


0.045 


3.0 x 10- 


3 


1 x 10" 


-3 


31 


0.450 


3.0 x 10" 


2 


6 x 10" 


-3 


25 


2.700 


1.8 x 10- 


1 


1 x 10" 


-2 


19 


4.500 


3.0 x 10- 


1 



Table 1. Summary of potential strengths a/U used in the numerical simulations. 

Also given the maximum odd particle number N < min (2/ y/a/U , 31 J that docs 

not lead to unit filling in the trap core at J/U = for a system of 61 sites, the 
potential drop between the center and the edges of the system Afi/U, and the 
maximum potential difference between neighboring sites aR/U. 



at J/U = docs not possess a core region with unit filling. The four different 
potential strengths a/U give rise to a variation of the local chemical potential from the 
center to the edge of A/J./U = fi /U - /j, r /U = R 2 (a/U)R 2 /2 = 0.045,0.45,2.7,4.5, 
respectively (this is the length of the system line introduced in the preceding section) . 
The smallest value is much smaller than the extent of the DW domain in the phase 
diagram (Fig. [IJ, the intermediate ones are comparable, and the largest one is much 
bigger. The potential difference between neighboring sites remains smaller than 
{a/U)R = 3 • 10~ 3 ,3 • 10~ 2 ,1.8 • 10 _1 ,3 • 10 -1 , respectively. Hence, even for the 
steepest potential at tunneling strength J/U > 1/2 i.e. before entering the DW regime, 
the particles are still delocalized over several sites. The numbers presented in this 
paragraph are summarized in tabic [I] 

The system is initialized in its ground state for J/U = 0.7, before J/U is ramped 
down linearly from 0.7 to within the time span T = tH/U. We choose values 
between r = 1 (intermediate) and r = 10 (moderately large, even larger times would 
be desirable but are numerically costly). 

After the ramp, we compute the distance to the target state of a perfect DW 
crystal with exactly one particle on every other site in the central region of the trap. 
The first measure we consider for this purpose is the DW order parameter. However, 



we are not using Odw &s it is defined in Eq. (18), but instead the definition 



a £<,ygAf {ntne) (-1)' " , . 

CdW = "2 (/Uj 

which is well defined also for a region M' of finite extent L' only. For a perfect zig-zag 
(DW) structure 3/4 of the terms in the nominator is and the rest is 1. The sum in 
the denominator is L'/2 for a perfect zig-zag structure so the ratio becomes exactly 1 
in this case. In order to exclude edge effects and to be able to compare scenarios with 
different particle numbers, we compute Odw based on a region M' C M containing the 
central 31 sites (L' « L/2). In an experiment the density-density correlations (n^n^) 
entering this parameter can be extracted by site resolved measurements |61[ I62j . 

As a second measure we use the nearest-neighbor fidelity -Fnn- We compute the 
two-site reduced density matrix for each pair of neighboring sites I and I'. For the 
time evolved state \tp) this 4x4 matrix is defined by 

$ =TreeM\ { e,e> } mm- (21) 

The two-site reduced density matrix is sufficient to calculate all two-site observables, 
and therefore characterizes spatially local properties of the state. We can now calculate 
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the symmetrized overlap ^(PnniPdw) between the two-site density matrices of the 
time evolved state, ffl , and those computed for the target state with perfect DW 
order, Pdw- The symmetrized overlap between two density matrices reads 

d(p A , PB ) = Tr (^j y/pX P By/p^j (22) 

which is 1 if and only if pa = Pb [63] . The nearest-neighbor fidelity is then defined as 
the average of these overlaps over all neighboring sites in M': 

Jnn = L , _ r ^ ^(PnnjPdw)- ( 23 ) 

We plot the main results of our simulation in Fig. [2j Both measures the DW 
order parameter Odw and the nearest-neighbor fidelity -Fnn give the same qualitative 
picture. We find the best result (the largest degree of adiabaticity) for the system 
at "half filling" with N = 31 particles in combination with the shallowest parabolic 
potential. However, the degree of adiabaticity that has been achieved for generic 
particle numbers below half filling (N < 29) is almost as large as for half filling, and 
it can be increased further by tuning the potential depth a/U continuously to its 
optimal value for every particle number (instead of using only four different values 
of a/U as we do here). So we prefer not to emphasize the better results for half 
filling. We observe very clearly that the optimal trapping strength depends sensitively 
on the particle number; the optimal depth a/U increases when the particle number 
N is lowered. This effect tells us that a finite parabolic inhomogeneity generally 
does assist the adiabatic preparation of the DW crystal. As expected, the degree of 
adiabaticity increases with r; the tendency of the curves suggest that the results can 
still be considerably improved by using ramp times larger than r = 10. 

In order to get further insight, in Fig. [3] we report the time evolution of the system 
with N = 29 and a/U — 1CP 3 during the ramp with r = 10. Panel (a) and (b) show 
the single-particle correlations (StS ) and the density-density correlations (fieno) both 
before the ramp (solid lines) and in the middle of the ramp (dashed lines). As expected, 
we can observe that with time the single-particle correlations (the off-diagonal order) 
decrease whereas the density-density correlations of the DW type are increased. This 
behavior is also reflected in the fact that the DW order parameter Odw as well as the 
nearest-neighbor fidelity Fnn grow with time [panel (c) and (d)]. 

A more subtle effect is visible in the density-density correlations shown in 
Fig. |3jb). Already for the initial state (the ground state at J/U = 0.7, see solid 
line) it posseses traces of a DW-type zig-zag pattern. Superimposed to this pattern 
one can observe a modulation on a larger length scale (comparable to the system size) 
having nodes roughly at £ = ±16. Having a closer look reveals that at these nodes 
the zig-zag correlations have a kink. Now, it is very difficult for the system to get 
rid of the kinks during the ramp as it would be required for a perfectly adiabatic 
evolution (since eventually at J/U = the ground state is a perfect DW). Therefore 
during the ramp the kinks remain, i.e. they are converted into defects. This becomes 
evident from Fig. |4ja) where the density distribution of the time-evolved state after 
the ramp is plotted (the other subfigures of Fig. [4] display more information on the 
final state). The presence of the kinks explains also the significant drop of the DW 
order parameter Odw when computed not only for 31 central sites but rather for the 
whole system [Fig. [3jc)] . 
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Figure 2. Degree of adiabaticity during a ramp from the supcrfiuid to the 
density wave (DW) regime for a system of N particles on a lattice of L = 61 
sites in the presence of a parabolic potential of strength a/U = 10~ 4 (blue dotted 
lines), 10 -3 (dashed red lines), 6 • 10 — 3 (dash-dotted green lines up to N = 25), 
10 -2 (dash-dotted purple lines up to N = 19). Starting from the ground state 
at a tunneling parameter of J/U = 0.7, the time evolution is simulated while 
J/U is linearly ramped down to zero within a time span tH/U with r = 3 
(crosses), 5 (circles), 7 (diamonds), 9 (squares), 10 (triangles). For the final 
state we plot (a) the normalized DW order parameter Odw an< 3 (t>) the nearest- 
neighbor fidelity Fnjm with respect to the DW ground state (b). Both quantities 
are computed for the central region of 31 sites and approach unity in the limit 
of a perfectly adiabatic dynamics. The best results are found for the largest 
particle number N = 31 (corresponding to "half filling" ) in combination with the 
shallowest parabolic potential. In contrast, for the lower particle numbers N < 29 
the presence of steeper potentials is always found to be favorable. As expected, 
the degree of adiabaticity increases with r; the tendency of the curves suggest 
that the results can be improved further by using ramp times larger than r = 10. 
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Figure 3. Sample time evolution of a system of N = 29 particle on L = 61 sites 
in the presence of the steep parabolic potential of strength a/U = 10 — 3 . The 
dimensionless tunneling parameter J/U is linearly ramped from 0.7 to over a 
duration of r = 10 (in units of h/U). (a) Single-particle density matrix {b\bj) for 
j = at times t = (solid line) and t = r/2 = 5 (dashed line), (b) Density- 
density correlations (hihj) for j = at t = and t = r/2 = 5. (c) Time evolution 
of the order parameter Odw f° r the entire system and for a 30 site cut-out of the 
central part of the trap, (d) Average of the time evolution of the nearest-neighbor 
density matrix fidelity Fnn for the central part of the trap. 
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Figure 4. Final state of a system with N = 29 particles on L = 61 after a 
ramp of J/U from 0.7 to with ramp duration r = 10 and potential strength 
a/E/ = 10~ 3 . (a) Site occupation (hi) (b) Density-density correlations {niUj). 
(c) Single-site particle number variance m = (n?) — ("i} 2 (d) Diagonal momentum 
correlation function (o^ol; ay Ofe) as it can be extracted from the shot noise 
correlations of time-of-flight absorption images. The satellite peaks at 0.5 are 
a signature of the density-wave order. 



For other particle numbers and potential strengths we find similar results. Namely 
the initial ground state at J/U = 0.7 possesses already a small DW-type modulation of 
the site occupation, typically contaminated with superimposed large scale modulations 
and/or a few kinks. These initial density correlations, including the kinks, are 
amplified when J/U is lowered within a time of rh/U. This behavior can be inferred 
from Fig. [5] that shows the initial and the final density distribution for several particle 
numbers and trap depths. The best (most adiabatic) results are found when there are 
no kinks or if the kinks lie outside the central 31-site region that we use to measure 
the DW order. 
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N = 19 N = 19 N = 23 N = 25 N = 31 JV = 31 

a/C/ = 1 X 10~ 2 a/(7 = 6 X 10~ 3 a/(7 = 6 X 10~ 3 a/(7 = 6 X 10 -3 a/U = 1 X 10~ 3 a/U = 1 X 10~ 4 




Figure 5. Density distribution before and after a ramp of duration rh/U 
with t — 10 for different particle number and trap depths. The low degree of 
adiabaticity for N = 19 with a/U = 6 ■ 10 — 3 can be ascribed to the structure of 
the initial state. It possesses already weak DW correlations that are contaminated 
with defects. The system is not able to get rid of these defects during the ramp. 
This behavior is found for trap depth (particle numbers) that are weaker (larger) 
than optimal. The low degree of adiabaticity for both N = 25 with a/U = 6- 10 — 3 
and 31 with a/U = 1T0~ 3 is related to density modulations on larger scales. This 
behavior is found for trap depth (particle numbers) that are stronger (smaller) 
than optimal. 

We find that the results are spoiled by kinks typically when the trapping potential 
is (according to Fig. [2| shallower than optimal for a given particle number (or the 
particle number is larger than optimal for a given potential strength). A typical 
example for kinks spoiling an adiabatic time evolution is found for N = 19 with 
a/U = 6 ■ 10~ 3 (Fig. [5]). If, in turn, the trap is to steep for a given particle number (or 
the particle number too low for a given trap depth) , we observe superimposed density 
modulations on larger scales in the initial state (not necessarily in combination with 
kinks). These modulations are still found in the time evolved state after the ramp, 
so they lower the degree of adiabaticity. In Fig. [5] this behavior is visible for TV = 25 
with a/U = 6 • 1(T 3 and N = 31 with a/U = 1 • 1(T 3 . 

The superimposed large-scale modulation of the DW correlations as well as the 
kinks that are present in the initial state and hamper an adiabatic time-evolution 
during the ramp cannot be explained within the simple picture of the LDA. They 
originate from the trap and the finite extent of the system. This suggests that the 
picture that in the previous section was drawn on the basis of the LDA (augmented 
by the assumption of smooth crossovers at phase boundaries) does not yet apply 
completely for the experimentally relevant system sizes of 50-100 sites only. 

One might speculate that kinks and density modulations are an artifact of the 
open boundary condition. However, as becomes apparent in Fig. [5| we also find 
kinks in the initial states for the steep trapping potentials with a/U — 6 • 10 -3 and 
a/U — 1-10~ 2 for which the initial state has practically no occupation at the outermost 
sites such that the boundary conditions do not matter. Nevertheless, for the shallower 
trapping potentials, the initial state does depend on the artificial open boundary 
conditions and the finite size effects that we observe here can be modified when the 
edge of the mixed Mott insulator domain M is not approximated by open boundary 
conditions. A more realistic model that captures also the shell surrounding the mixed 
Mott-insulator region M is given by Eq. ([3]). In Fig. kS] we plot the occupation numbers 




Figure 6. Site atom numbers for the two-species model in the regime used to 
simulate the XXZ-model with N a = 39 and JVj, = 28 particles. The parameters 
are U a b = 1.0, ct a /b ~ 9.1 X 10 _4 ±6.0 X 10~ 7 . The two-species system has 81 sites 
and effective tunneling corresponding to (a) J/U = 0.5 (J a = Jb = U a b/&) and 
(b) J/U = 0.12 (J a = U aa /S, Jb = J a /8). The trap inhomogeneity correspond to 
approximately a = 10 — 6 . 



of a and b particles for the ground state of this two-species model at J/U = 0.5 
and J/U = 0.1 (the other parameters are specified in the caption). For the small 
tunneling the ground state features defect- free DW/antiferromagnetic order in the 
central Mott region. The small antiferromagnetic correlations found for the larger 
tunneling parameter are, however, contaminated with defects in a similar way as 
observed before when the Mott region M with open boundary conditions was treated. 
These defects can spoil an adiabatic parameter variation towards the antiferromagnetic 
state plotted in Fig. fBVb). 

In an experiment the DW order can be observed either in situ using high- 
resolution imaging techniques |61j , or via time-of- flight noise correlation measurements 
[64l 165"! 166] probing the two-particle momentum correlation function plotted in 
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Fig. Qd). In the latter, the signature of the DW order is given by the two satellite 
peaks. The fact that this feature is rather small is a consequence of the fact that in 
our simulations the ramp times are still not large enough. For larger ramp times (the 
simulation of which is numerically costly) the degree of adiabaticity is still expected 
to improve considerably. 

6. Conclusion and outlook 

We have pointed out that in a mixed Mott insulator of two bosonic atom species 
in an optical lattice a parabolic inhomogeneity can be created and widely tuned 
(between zero and large finite values) by introducing a finite potential difference for 
both species. And we proposed to use this control knob to investigate the role of such 
an inhomogeneity on the adiabatic preparation of an antiferromagnetic state (with 
a staggered DW pattern for each of the species). Numerical simulations of the time 
evolution of a model describes such a system in one dimension lead to the conclusion 
that a finite inhomogeneity generally assists the adiabatic preparation of the bosonic 
antiferromagnet. The optimal strength of the parabolic inhomogeneity depends in 
a sensitive way on the imbalance between the particle numbers of the two species; 
the larger the imbalance the larger the optimal strength of the inhomogeneity. We 
find that for a realistic system size (a mixed Mott insulator stretched over 60 sites) 
finite size effects that cannot be explained within the local density approximation are 
significant. Namely the mechanism leading to deviations from adiabaticity is related 
to the presence of precursing DW modulations already outside the antiferromagnetic 
parameter regime that - as a consequence of the finite system size - comprise 
imperfections like kinks. When the system is ramped into the antiferromagnetic regime 
these modulations, including the imperfections, are amplified. 

We believe that the experimental implementation of tunable parabolic potential 
as we propose it here, can be a valuable tool for finding a good protocol for the 
preparation of antiferromagnetic order. The concept generalizes also to two and three 
spatial dimensions, a situation which can be addressed easily in an experiment. A 
theoretical study of the higher-dimensional case could be carried out on a qualitative 
level using Gutzwiller mean-field theory. 

Another relevant question would be whether such a controllable parabolic 
potential can be useful also for the preparation of antiferromagnetic order in a Mott 
insulator of fermionic atoms [671 168) . 
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Appendix A. Bethe-Ansatz solution of the homogeneous problem 



The eigenstates of the Hamiltonian (17) with a homogeneous magnetic field he = h is 
found in Refs. [55J [SHI [57] using the Bethe Ansatz. 
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Let y be the total z-magnetization density and A = — U/2J < characterizes the 
interaction strength. Then, A = — 1 correspond to the Heisenberg anti-ferromagnet 
as mentioned above. 



Instead of finding the groundstate of (17) for a specific magnetic field we will 
instead find the ground state in each total spin-z subspace, or equivalently, for fixed 
particle number. To find the ground state at a given magnetic field h, or chemical 
potential fi = h we should then minimize u(A, y) — \iy, where u = E/L is the energy 
density of the ground state, with respect to y, giving /i = du/dy(A,y). 

In [57] the phase boundary between the SF and DW phases in Fig. [I] is given 
analytically by 

-^=2sinh(A) V \ / . , (A.l) 
2 J 1 ' ^ cosh(nA) V ' 

n— — oo v * 

where cosh(A) = -A = U/2J. 

The ground state energy density u = E/L and magnetization < y < 1 can be 
found by solving the integral equations [561 Eq. (7a-c)] 

R(a) =^r- t K(a,P)R(B)dp (A.2) 
da J_ b 

\)da (A.3) 

2 cL R ^ta da ^ 

where b G [0, tt] for A < — 1 and b £ [0, oo) for — 1 < A < 1 and the functions K and 
dp /da is given in [56L Table II]. 

For fixed 6, the first equation is a Fredholm integral equation of the second kind, 
and it can be solved using the Nystrom method [69] • For fixed b we can therefore solve 
for R and, using this solution, calculate the magnetization y using the second equation 
and the energy density u using the third. The magnetic field is fj,/U — du/dy(A,y) 
which can also be found by calculating y for b ± db for a small db numerically and 
approximating the derivative using the finite difference. 

One can show that the function R is positive, which implies the y is a one-to- 
one function of b. In order to determine the phase diagram in Fig. [T] we calculate 
y and /i = de/dy as a function of A and b. This results in a set of (non- uniformly 
distributed) points (A, b, y, fJ-/U), where we can now plot y as a function of (A, fi/U), 
ov,(J/U,n/U). 



(i-y) = 


[ R{ 




J-b 


(A,y) 


A 


2J 


~ ~2 ~ 
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